Introduction

We analysed 1434 flights of 50 individuals of blue tit, great tit, and marsh tit as baseline data We found significant difference in initial velocity and in-flight acceleration between species . Great tits showed higher initial velocity and lower in-flight acceleration than blue tits. Marsh tits did not differ in flight performance with other species. Weather conditions did not show significant effect to flight performance.

I. Data

libraries

library(dplyr)
library(stringr)
library(readr)
library(fuzzyjoin)
library(tidyverse)
library(ggplot2)
library(lubridate)
library(ggpubr)
library(rstatix)
library(coin)
library(scales)
library(gridExtra)
library(ggthemes)
library(ggbreak)
library(lme4)
library(reshape2)
library(olsrr)
library(gvlma)
library(PerformanceAnalytics)
library(usethis)
library(devtools)
library(redres)
library(EnvStats)
library(lmtest)
library(lmerTest)
library(AICcmodavg)
library(emmeans)
library(MuMIn)
library(sjPlot)

dataset

#in dept laptop
#setwd("C:\\Users\\jesu4036\\OneDrive - Nexus365\\Documents\\R\\WythamWoods")
#load all dates data 
#all_dates3 <- read.delim(file = ".\\all_dates3.txt", sep = ",")

#in the kyubird-laptop
all_dates3 <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\all_dates3.txt", sep = ",") %>% distinct()

all_dates3$datetime.y <- as.POSIXct(all_dates3$datetime.y, format = "%Y-%m-%d %H:%M:%S")
colnames(all_dates3)

##  [1] "Ring"        "RFID"        "datetime.y"  "temperature" "firstbreak1"
##  [6] "firstbreak2" "firstbreak3" "lastbreak1"  "lastbreak2"  "lastbreak3" 
## [11] "t1"          "t2"          "v1"          "v2"          "avg.v"      
## [16] "acc"         "section2"    "section3"    "section4"    "length1"    
## [21] "length2"     "tunnel"      "species"     "age"         "Spec"       
## [26] "Sex"         "Date"        "Place"       "Site"        "Wing"       
## [31] "Wt"          "Init"        "Tarsus"

all_dates3 <- all_dates3 %>% 
  select(-Spec)
 

#load climate data 
#in dept laptop 
#climate <- read.delim(file = ".\\climate_filled.csv", sep = ',')

#in the kyubird laptop
climate <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\climate_filled.csv", sep = ',')

climate$datetime <- as.POSIXct(climate$datetime, format = "%Y-%m-%d %H:%M:%S")

#in the dept laptop 
#caught.hour.day <- read.delim(file = ".\\caught.hour.day.csv", sep = ",")

#in the kyubird laptop
caught.hour.day <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\caught.hour.day.csv", sep = ',')

adding day and hour and climate

all_dates3$day <- day(all_dates3$datetime.y)
all_dates3$hour <- hour(all_dates3$datetime.y)

all_dates3$roundate <- round_date(all_dates3$datetime.y, "10 mins")
all_dates3.climate <- merge(x = all_dates3, y = climate, by.x = "roundate", by.y = "datetime", all.x = TRUE)

flight.data <- all_dates3.climate %>%
  filter(day %in% c(15:24))

nat.var <- all_dates3.climate %>%
  filter(day %in% c(15:19))

capt.eff <- all_dates3.climate %>%
  filter(day %in% c(20:24))

write.csv(nat.var, “C:\Users\kmh\Desktop\DATA\nat.var.csv”)

data synopsis

raw continuous data

dim(flight.data) #total number of flight

## [1] 3451   38

table(flight.data$day)

## 
##  15  16  17  18  19  20  21  22  23  24 
## 241 231  54 528 380 485 616 318 504  94

no.of.individuals <- rowSums(table(flight.data$species, flight.data$Ring)!=0) 
no.of.flights <- table(flight.data$species)

raw_sum <- rbind(no.of.individuals, no.of.flights)

raw_sum

##                   bluti greti marti
## no.of.individuals    40    16     5
## no.of.flights      2617   513   321

natural variation data

dim(nat.var) #total number of flight

## [1] 1434   38

table(nat.var$day)

## 
##  15  16  17  18  19 
## 241 231  54 528 380

no.of.individuals <- rowSums(table(nat.var$species, nat.var$Ring)!=0) 
no.of.flights <- table(nat.var$species)

nat.var_sum <- rbind(no.of.individuals, no.of.flights)

nat.var_sum

##                   bluti greti marti
## no.of.individuals    33    12     4
## no.of.flights      1172   148   114

breeding <- cbind(nat.var$Ring, nat.var$RFID, nat.var$species, nat.var$age, nat.var$Sex) %>% as.data.frame() %>%
  distinct()
colnames(breeding) <- c("Ring", "RFID", "species", "age", "Sex")

write.csv(breeding, “\breeding.csv”,row.names = FALSE )

For each species, M vs. F, young vs. adult, then between species

#find the ring number of those without Sex in Great Tit & Blue Tit 
for( i in 1:nrow(nat.var)){
  
  if(nat.var$Ring[i] == "ATR8387"){
    
    nat.var$Sex[i] <- "F"
  }
  
  else if(nat.var$Ring[i] == "ATR8383"){
    
    nat.var$Sex[i] <- "F"
  }
  else if(nat.var$Ring[i] == "PF72506"){
    
    nat.var$Sex[i] <- "M"
  }
  else{}
  
}

#Sex unknown for ATR8380 and ATR5368 both in morphometrics and brood data
nosex <- nat.var %>%
  filter(species != "marti") %>%
  filter(Sex == "") %>%
  select(Ring) %>%
  distinct()

no.of individuals by species, age, sex in natural variation

numbers <- nat.var %>%
  count(species, age, Sex, Ring)


individual.age <- numbers %>%
  count(species, age)

individual.sex <- numbers %>%
  count(species, Sex)

numbers.age <- numbers %>%
   group_by(species, age) %>%
   summarise(mean = mean(n), sd = sd(n)) %>%
             na.omit()

## `summarise()` has grouped output by 'species'. You can override using the
## `.groups` argument.

numbers.sex <- numbers %>%
  group_by(species, Sex) %>%
  summarise(mean = mean(n), sd = sd(n)) %>%
            na.omit()

## `summarise()` has grouped output by 'species'. You can override using the
## `.groups` argument.

numbers.age <- numbers.age %>%
  add_column(individual.age$n)

numbers.sex <- numbers.sex %>%
  add_column(individual.sex$n)

numbers.tot <- numbers %>%
  group_by(species) %>%
  summarise(mean = mean(n), sd = sd(n)) %>%
            na.omit()

numbers.age

## # A tibble: 5 × 5
## # Groups:   species [3]
##   species   age  mean    sd `individual.age$n`
##   <chr>   <int> <dbl> <dbl>              <int>
## 1 bluti       5  33.6  24.6                 22
## 2 bluti       6  36    20.0                 12
## 3 greti       5  11.4  11.5                 10
## 4 greti       6  17    22.6                  2
## 5 marti       5  28.5  11.0                  4

numbers.sex

## # A tibble: 6 × 5
## # Groups:   species [3]
##   species Sex    mean    sd `individual.sex$n`
##   <chr>   <chr> <dbl> <dbl>              <int>
## 1 bluti   ""     54.5  17.7                  2
## 2 bluti   "F"    33.7  18.9                 15
## 3 bluti   "M"    32.8  26.2                 17
## 4 greti   "F"    15.8  17.4                  4
## 5 greti   "M"    10.6  10.5                  8
## 6 marti   ""     28.5  11.0                  4

numbers.tot

## # A tibble: 3 × 3
##   species  mean    sd
##   <chr>   <dbl> <dbl>
## 1 bluti    34.5  22.8
## 2 greti    12.3  12.6
## 3 marti    28.5  11.0

t-tests on visiting rate

bluti <- numbers %>%
  filter(species == "bluti")

greti <- numbers %>%
  filter(species == "greti")

marti <- numbers %>%
  filter(species == "marti")



bluti.F <- numbers %>%
   filter(species == "bluti" & Sex == "F")

bluti.M <- numbers %>%
   filter(species == "bluti" & Sex == "M")

bluti.5 <- numbers %>%
   filter(species == "bluti" & age == 5)

bluti.6 <- numbers %>%
   filter(species == "bluti" & age == 6)

greti.F <- numbers %>%
  filter(species == "greti" & Sex == "F")

greti.M <- numbers %>%
  filter(species == "greti" & Sex == "M")

greti.5 <- numbers %>%
  filter(species == "greti" & age == 5)

greti.6 <- numbers %>%
  filter(species == "greti" & age == 6)

t.test(bluti.F$n, bluti.M$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti.F$n and bluti.M$n
## t = 0.10527, df = 28.948, p-value = 0.9169
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -15.53856  17.22483
## sample estimates:
## mean of x mean of y 
##  33.66667  32.82353

t.test(bluti.5$n, bluti.6$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti.5$n and bluti.6$n
## t = -0.30328, df = 27.045, p-value = 0.764
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -18.35356  13.62629
## sample estimates:
## mean of x mean of y 
##  33.63636  36.00000

t.test(greti.F$n, greti.M$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  greti.F$n and greti.M$n
## t = 0.5423, df = 4.1323, p-value = 0.6155
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -20.7855  31.0355
## sample estimates:
## mean of x mean of y 
##    15.750    10.625

t.test(greti.5$n, greti.6$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  greti.5$n and greti.6$n
## t = -0.34135, df = 1.1049, p-value = 0.7863
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -172.6219  161.4219
## sample estimates:
## mean of x mean of y 
##      11.4      17.0

t.test(bluti$n, greti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti$n and greti$n
## t = 4.146, df = 35.353, p-value = 0.0002014
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  11.30160 32.97291
## sample estimates:
## mean of x mean of y 
##  34.47059  12.33333

t.test(bluti$n, marti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti$n and marti$n
## t = 0.88326, df = 6.6185, p-value = 0.408
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -10.20223  22.14341
## sample estimates:
## mean of x mean of y 
##  34.47059  28.50000

t.test(greti$n, marti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  greti$n and marti$n
## t = -2.447, df = 5.8748, p-value = 0.05085
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -32.41668686   0.08335352
## sample estimates:
## mean of x mean of y 
##  12.33333  28.50000

headwind, air temperature, humidity

#overall climate data 
#find how many weather values are missing
#missing data that do not make to the final dataset. 
View(flight.data[!complete.cases(flight.data),])
#mean weather 

mean(flight.data$air.temp, na.rm = TRUE)

## [1] 9.096442

mean(flight.data$headwind, na.rm = TRUE) #relative to the birds flying out of the tunnel

## [1] -1.428632

mean(flight.data$humidity, na.rm = TRUE)

## [1] 80.46904

#dataset for baseline variation 
climate.nat.var <- nat.var %>%
  filter(day %in% c(15:19)) %>%
  select(headwind, air.temp, humidity)
  
#not using caught.hour.day which is a reduced dataset to just answer the question if the climate conditions were different at nat.var vs. capt.eff
#otherwise everything is significant probably the sample size is so vastly different
climate.capt.eff <- caught.hour.day %>%
  select(headwind, air.temp, humidity)

#mean of nat.var. # 2 for columns
apply(climate.nat.var, 2,mean)

##  headwind  air.temp  humidity 
## -3.169797 10.976903 78.366511

#sd of nat.var
apply(climate.nat.var, 2, sd)

## headwind air.temp humidity 
## 1.040167 1.553028 9.631564

#mean of capt.eff
apply(climate.capt.eff, 2, mean, na.rm = TRUE)

##  headwind  air.temp  humidity 
##  2.300897  6.819277 72.716639

#mean of capt.eff
apply(climate.capt.eff, 2, sd, na.rm = TRUE)

##  headwind  air.temp  humidity 
## 2.1071941 0.9258642 5.7299068

climate.all <- nat.var %>%
  filter(day %in% c(15:24)) %>%
  select(headwind, air.temp, humidity)

apply(climate.all, 2, mean)

##  headwind  air.temp  humidity 
## -3.169797 10.976903 78.366511

apply(climate.all, 2, sd)

## headwind air.temp humidity 
## 1.040167 1.553028 9.631564

caught.byday <- caught.hour.day %>%
  group_by(cgt_day) %>%
  summarise(mean.hdw = mean(headwind, na.rm = FALSE),
            mean.art = mean(air.temp, na.rm = FALSE),
            mean.hmd = mean(humidity, na.rm = FALSE),
            sd.hdw = sd(headwind, na.rm = FALSE),
            sd.art = sd(air.temp, na.rm = FALSE),
            sd.hmd = sd(humidity, na.rm = FALSE),
            max.hdw = max(headwind, na.rm = FALSE),
            min.hdw = min(headwind, na.rm = FALSE)
            ) %>%
            na.omit()


caught.byday

## # A tibble: 2 × 9
##   cgt_day mean.hdw mean.art mean.hmd sd.hdw sd.art sd.hmd max.hdw min.hdw
##     <int>    <dbl>    <dbl>    <dbl>  <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
## 1       0    0.156     7.58     75.8  1.53   0.774   5.46    4.83   -1.05
## 2       1    3.71      6.32     70.7  0.894  0.633   4.98    5.47    1.99

day0 <- caught.hour.day %>%
          filter(cgt_day == 0)

day1 <- caught.hour.day %>%
          filter(cgt_day == 1)

t.test(day0$headwind, day1$headwind, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  day0$headwind and day1$headwind
## t = -18.211, df = 107.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.942435 -3.168404
## sample estimates:
## mean of x mean of y 
## 0.1563588 3.7117780

t.test(day0$air.temp, day1$air.temp, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  day0$air.temp and day1$air.temp
## t = 11.723, df = 136.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.045330 1.469544
## sample estimates:
## mean of x mean of y 
##  7.577731  6.320294

t.test(day0$humidity, day1$humidity, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  day0$humidity and day1$humidity
## t = 6.4736, df = 148.01, p-value = 1.323e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.527083 6.626559
## sample estimates:
## mean of x mean of y 
##  75.77885  70.70203

t.test(climate.nat.var$headwind, climate.capt.eff$headwind, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  climate.nat.var$headwind and climate.capt.eff$headwind
## t = -35.132, df = 200.24, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.777751 -5.163637
## sample estimates:
## mean of x mean of y 
## -3.169797  2.300897

t.test(climate.nat.var$air.temp, climate.capt.eff$air.temp, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  climate.nat.var$air.temp and climate.capt.eff$air.temp
## t = 52.727, df = 347.02, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  4.002540 4.312713
## sample estimates:
## mean of x mean of y 
## 10.976903  6.819277

t.test(climate.nat.var$humidity, climate.capt.eff$humidity, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  climate.nat.var$humidity and climate.capt.eff$humidity
## t = 11.571, df = 347.77, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  4.689545 6.610200
## sample estimates:
## mean of x mean of y 
##  78.36651  72.71664

duration difference

raw.flight.metrics <- nat.var %>%
  group_by(species) %>%
  summarise(mean.v1 = mean(v1, na.rm = FALSE),
            mean.v2 = mean(v2, na.rm = FALSE),
            mean.acc = mean(acc, na.rm = FALSE),
            mean.avg.v = mean(avg.v, na.rm = FALSE),
            mean.t1 = mean(t1, na.rm = FALSE),
            mean.t2 = mean(t2, na.rm = FALSE),
            mean.t3 = mean((t1+t2), na.rm = FALSE),
            sd.v1 = sd(v1, na.rm = FALSE),
            sd.v2 = sd(v2, na.rm = FALSE),
            sd.acc = sd(acc, na.rm = FALSE),
            sd.avg.v = sd(avg.v, na.rm = FALSE),
            sd.t1 = sd(t1, na.rm = FALSE),
            sd.t2 = sd(t2, na.rm = FALSE),
            sd.t3 = sd((t1+t2), na.rm = FALSE)
            ) %>%
            na.omit()

raw.flight.metrics

## # A tibble: 3 × 15
##   species mean.v1 mean.v2 mean.acc mean.avg.v mean.t1 mean.t2 mean.t3 sd.v1
##   <chr>     <dbl>   <dbl>    <dbl>      <dbl>   <dbl>   <dbl>   <dbl> <dbl>
## 1 bluti      2.53    2.94   3.89         2.71   107.     105.    212. 0.451
## 2 greti      2.84    2.87  -0.0160       2.83    97.4    108.    206. 0.624
## 3 marti      2.62    2.85   1.96         2.71   106.     108.    214. 0.511
## # ℹ 6 more variables: sd.v2 <dbl>, sd.acc <dbl>, sd.avg.v <dbl>, sd.t1 <dbl>,
## #   sd.t2 <dbl>, sd.t3 <dbl>

bluti.raw <- nat.var %>%
  filter(species == "bluti") %>%
  mutate(t3 = t1 + t2)

greti.raw <- nat.var %>%
  filter(species == "greti") %>%
  mutate(t3 = t1 + t2)





t.test(bluti.raw$t1, greti.raw$t1, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti.raw$t1 and greti.raw$t1
## t = 3.943, df = 192.45, p-value = 0.0001126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   4.857541 14.581072
## sample estimates:
## mean of x mean of y 
## 107.13823  97.41892

t.test(bluti.raw$t2, greti.raw$t2, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti.raw$t2 and greti.raw$t2
## t = -1.6142, df = 207.97, p-value = 0.108
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.0957960  0.7070403
## sample estimates:
## mean of x mean of y 
##  104.9881  108.1824

t.test(bluti.raw$t3, greti.raw$t3, paired = FALSE, var.equal = FALSE, conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  bluti.raw$t3 and greti.raw$t3
## t = 1.6545, df = 188.86, p-value = 0.09969
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.25462 14.30448
## sample estimates:
## mean of x mean of y 
##  212.1263  205.6014

write.csv(raw.flight.metrics, “\raw.flight.metrics.csv”,row.names = FALSE )

II. Graphs

A. Structural variation

1. individual variation

2. species variation

3. sex x species variation

4. age variation

B. climate factors

5.1 headwind

5.2 air temperature

5.3 humidity

C. diurnal variation

day effect not included in linear or linear mixed models because there was no difference among days until mist-netting happened.

III. Linear Mixed Model analysis

We chose linear mixed model to account for individual variation in flight speed. We chose headwind from climate factors by its relationship with flight and correlation with air temperature and humidity.

A. Model selection

1. checking correlations

##  [1] "roundate"    "Ring"        "RFID"        "datetime.y"  "temperature"
##  [6] "firstbreak1" "firstbreak2" "firstbreak3" "lastbreak1"  "lastbreak2" 
## [11] "lastbreak3"  "t1"          "t2"          "v1"          "v2"         
## [16] "avg.v"       "acc"         "section2"    "section3"    "section4"   
## [21] "length1"     "length2"     "tunnel"      "species"     "age"        
## [26] "Sex"         "Date"        "Place"       "Site"        "Wing"       
## [31] "Wt"          "Init"        "Tarsus"      "day"         "hour"       
## [36] "air.temp"    "headwind"    "humidity"

2. Proposed Mixed Linear Models.

with headwind

lmer.v1 = lmer(v1 ~ species  + Sex  + age + headwind + (1|Ring), data = nat.var)
lmer.v2 = lmer(v2 ~species  + Sex  + age + headwind + (1|Ring), data = nat.var)
lmer.avg.v = lmer(avg.v ~ species  + Sex  + age + headwind + (1|Ring), data = nat.var)
lmer.acc = lmer(acc ~ species  + Sex  + age + headwind + (1|Ring), data = nat.var)

with air temperature

lmer.v1 = lmer(v1 ~ species  + Sex  + age + air.temp + (1|Ring), data = nat.var)
lmer.v2 = lmer(v2 ~species  + Sex  + age + air.temp + (1|Ring), data = nat.var)
lmer.avg.v = lmer(avg.v ~ species  + Sex  + age + air.temp + (1|Ring), data = nat.var)
lmer.acc = lmer(acc ~ species  + Sex  + age + air.temp + (1|Ring), data = nat.var)

3. preparing dataset without outliers

We then identified 15 outliers using Rosner’s Test(1983) and prepared a dataset excluding outliers to see the sensitivity of the models to the outliers. All the outliers were included in the dataset in the end as we could not confirm that these outliers were from measurement errors.

test.v1 <- rosnerTest(nat.var$v1, k=3)
test.v1$all.stats

test.v2 <- rosnerTest(nat.var$v2, k=10)
test.v2$all.stats

test.avg.v<- rosnerTest(nat.var$avg.v, k=13)
test.avg.v$all.stats

test.acc <- rosnerTest(nat.var$acc, k=3)
test.acc$all.stats

#collecting outlier row numbers
out.v1 <- test.v1$all.stats %>% 
  filter(Outlier == "TRUE") %>%
  select(Obs.Num)

out.v2 <- test.v2$all.stats %>% 
  filter(Outlier == "TRUE") %>%
  select(Obs.Num)

out.avg.v <- test.avg.v$all.stats %>% 
  filter(Outlier == "TRUE") %>%
  select(Obs.Num)

out.acc <- test.acc$all.stats %>% 
  filter(Outlier == "TRUE") %>%
  select(Obs.Num)

outliers

outliers <- rbind(out.v1, out.v2, out.avg.v, out.acc) %>% distinct()
outliers #15 outliers

##    Obs.Num
## 1     1014
## 2      243
## 3     1318
## 4      951
## 5     1088
## 6     1177
## 7       27
## 8     1333
## 9      351
## 10     861
## 11    1222
## 12    1257
## 13    1207
## 14    1315
## 15     195

without outliers

nat.var.out <- nat.var[-c(outliers$Obs.Num),]
dim(nat.var.out) #1419 observations

## [1] 1419   38

B. Linear mixed model assumptions

normality tests

# Computes the default residuals (raw conditional)
raw_cond.v1 <- compute_redres(lmer.v1)

## Loading required namespace: testthat

raw_cond.v2 <- compute_redres(lmer.v2)
raw_cond.avg.v <- compute_redres(lmer.avg.v)
raw_cond.acc <- compute_redres(lmer.acc)

initial velocity

shapiro.test(raw_cond.v1) #p <0.05

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v1
## W = 0.9694, p-value < 2.2e-16

end velocity

shapiro.test(raw_cond.v2) #p <0.05

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v2
## W = 0.89926, p-value < 2.2e-16

average velocity

shapiro.test(raw_cond.avg.v) #p <0.05

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.avg.v
## W = 0.89764, p-value < 2.2e-16

end velocity

shapiro.test(raw_cond.acc) #p <0.05

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.acc
## W = 0.99357, p-value = 7.024e-06

checking assumption graphs

initial velocity

plot(nat.var$headwind, nat.var$v1) #linearity

plot(nat.var$air.temp, nat.var$v1) #linearity

hist(raw_cond.v1) #normality in histogram

plot_redres(lmer.v1) #heteroskedasticity

plot_resqq(lmer.v1) #normality of residuals

plot_ranef(lmer.v1) #normality of random effects residuals

end velocity

plot(nat.var$headwind, nat.var$v2)

plot(nat.var$air.temp, nat.var$v2) #wider range

hist(raw_cond.v2) #normality in histogram

plot_redres(lmer.v2) #heteroskedasticity

plot_resqq(lmer.v2) #normality of residuals. not great..

plot_ranef(lmer.v2) #normality of random effects residuals

average velocity

plot(nat.var$headwind, nat.var$avg.v)

plot(nat.var$air.temp, nat.var$avg.v)

hist(raw_cond.avg.v) #normality in histogram

plot_redres(lmer.avg.v) #heteroskedasticity

plot_resqq(lmer.avg.v) #normality of residuals. #not great..

plot_ranef(lmer.avg.v) #normality of random effects residuals

acceleration

plot(nat.var$headwind, nat.var$acc)

plot(nat.var$air.temp, nat.var$acc) #better linearity

hist(raw_cond.acc) #normality in histogram

plot_redres(lmer.acc) #heteroskedasticity

plot_resqq(lmer.acc) #normality of residuals. great.

plot_ranef(lmer.acc) #normality of random effects residuals

C. Model summary

initial velocity

library(lmerTest)
Anova(lmer.v1, type = 3, test.statistic = "F") #Kenward-Rodger df

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v1
##                   F Df  Df.res    Pr(>F)    
## (Intercept) 29.7686  1   40.38 2.671e-06 ***
## species      4.7090  2   39.58   0.01464 *  
## Sex          0.9514  2   38.38   0.39512    
## age          0.5195  1   39.58   0.47527    
## air.temp     0.8094  1 1425.45   0.36845    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.v1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v1 ~ species + Sex + age + air.temp + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1764.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7165 -0.5902  0.0500  0.6787  3.4932 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.05143  0.2268  
##  Residual             0.18525  0.4304  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.808e+00  5.135e-01  3.940e+01   5.467 2.76e-06 ***
## speciesgreti  2.844e-01  9.757e-02  4.739e+01   2.914  0.00543 ** 
## speciesmarti  2.374e-01  2.101e-01  3.235e+01   1.130  0.26669    
## SexF          7.792e-02  1.775e-01  3.166e+01   0.439  0.66366    
## SexM          1.724e-01  1.762e-01  3.167e+01   0.978  0.33531    
## age          -6.299e-02  8.719e-02  3.862e+01  -0.722  0.47436    
## air.temp     -7.243e-03  8.038e-03  1.425e+03  -0.901  0.36768    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr spcsmr SexF   SexM   age   
## speciesgret -0.219                                   
## speciesmart -0.444  0.049                            
## SexF        -0.418 -0.057  0.763                     
## SexM        -0.299 -0.149  0.742  0.895              
## age         -0.931  0.239  0.204  0.125 -0.004       
## air.temp    -0.161 -0.023 -0.008  0.003 -0.007 -0.008

lsmeans(lmer.v1, pairwise~species)

## $lsmeans
##  species lsmean     SE   df lower.CL upper.CL
##  bluti     2.47 0.0634 32.8     2.34     2.59
##  greti     2.75 0.1061 44.1     2.54     2.96
##  marti     2.70 0.1746 34.2     2.35     3.06
## 
## Results are averaged over the levels of: Sex, age 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate     SE   df t.ratio p.value
##  bluti - greti  -0.2844 0.0978 48.6  -2.906  0.0149
##  bluti - marti  -0.2374 0.2101 33.2  -1.130  0.5028
##  greti - marti   0.0469 0.2273 35.2   0.207  0.9768
## 
## Results are averaged over the levels of: Sex, age 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

end velocity

Anova(lmer.v2, type = 3, test.statistic = "F") #nothing significant

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v2
##                   F Df  Df.res    Pr(>F)    
## (Intercept) 45.5217  1   41.03 3.741e-08 ***
## species      0.0849  2   40.16    0.9188    
## Sex          0.0149  2   39.80    0.9852    
## age          0.5283  1   40.47    0.4715    
## air.temp     0.0280  1 1427.00    0.8671    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.v2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v2 ~ species + Sex + age + air.temp + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1283
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.7896 -0.4106  0.0900  0.5407  3.6279 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.04533  0.2129  
##  Residual             0.13154  0.3627  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   3.202e+00  4.736e-01  3.859e+01   6.760 4.83e-08 ***
## speciesgreti -2.248e-02  8.931e-02  4.566e+01  -0.252    0.802    
## speciesmarti -6.624e-02  1.955e-01  3.201e+01  -0.339    0.737    
## SexF          1.907e-03  1.654e-01  3.152e+01   0.012    0.991    
## SexM          1.414e-02  1.641e-01  3.152e+01   0.086    0.932    
## age          -5.862e-02  8.048e-02  3.806e+01  -0.728    0.471    
## air.temp      1.139e-03  6.792e-03  1.427e+03   0.168    0.867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr spcsmr SexF   SexM   age   
## speciesgret -0.221                                   
## speciesmart -0.447  0.049                            
## SexF        -0.420 -0.058  0.765                     
## SexM        -0.307 -0.148  0.746  0.899              
## age         -0.932  0.240  0.203  0.123  0.000       
## air.temp    -0.148 -0.022 -0.007  0.003 -0.006 -0.008

average velocity

Anova(lmer.avg.v, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: avg.v
##                   F Df  Df.res    Pr(>F)    
## (Intercept) 42.9523  1   41.16 6.934e-08 ***
## species      1.1445  2   40.28    0.3285    
## Sex          0.3973  2   40.09    0.6748    
## age          0.7669  1   40.64    0.3863    
## air.temp     0.1784  1 1426.98    0.6728    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.avg.v)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.v ~ species + Sex + age + air.temp + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1137.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5561 -0.4192  0.0955  0.6236  3.8317 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.04272  0.2067  
##  Residual             0.11863  0.3444  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   3.009e+00  4.582e-01  4.011e+01   6.566 7.48e-08 ***
## speciesgreti  1.253e-01  8.627e-02  4.729e+01   1.453    0.153    
## speciesmarti  9.603e-02  1.895e-01  3.334e+01   0.507    0.616    
## SexF          4.908e-02  1.603e-01  3.287e+01   0.306    0.761    
## SexM          1.023e-01  1.591e-01  3.287e+01   0.643    0.525    
## age          -6.834e-02  7.787e-02  3.960e+01  -0.878    0.385    
## air.temp     -2.729e-03  6.454e-03  1.427e+03  -0.423    0.672    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr spcsmr SexF   SexM   age   
## speciesgret -0.221                                   
## speciesmart -0.448  0.049                            
## SexF        -0.421 -0.058  0.766                     
## SexM        -0.308 -0.148  0.747  0.900              
## age         -0.933  0.240  0.203  0.122  0.000       
## air.temp    -0.145 -0.022 -0.007  0.003 -0.006 -0.008

acceleration

Anova(lmer.acc, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: acc
##                   F Df  Df.res    Pr(>F)    
## (Intercept)  3.1432  1   37.92   0.08428 .  
## species     13.1702  2   37.23 4.728e-05 ***
## Sex          0.4682  2   33.21   0.63018    
## age          0.6599  1   35.83   0.42195    
## air.temp     0.8469  1 1374.24   0.35758    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.acc)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: acc ~ species + Sex + age + air.temp + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 8283.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3928 -0.6556 -0.0021  0.6449  3.3918 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept)  2.119   1.456   
##  Residual             18.195   4.266   
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     6.44546    3.62535   33.99118   1.778   0.0844 .  
## speciesgreti   -3.43432    0.71333   43.71998  -4.814  1.8e-05 ***
## speciesmarti   -2.93246    1.42597   26.64274  -2.056   0.0497 *  
## SexF           -0.64160    1.19660   25.20720  -0.536   0.5965    
## SexM           -1.01659    1.18782   25.19098  -0.856   0.4002    
## age            -0.49848    0.61173   32.11682  -0.815   0.4211    
## air.temp        0.07240    0.07845 1368.22941   0.923   0.3562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr spcsmr SexF   SexM   age   
## speciesgret -0.201                                   
## speciesmart -0.428  0.047                            
## SexF        -0.406 -0.056  0.751                     
## SexM        -0.266 -0.149  0.725  0.880              
## age         -0.923  0.223  0.208  0.132 -0.020       
## air.temp    -0.224 -0.024 -0.010  0.000 -0.011 -0.011

lsmeans(lmer.acc, pairwise~ species)

## $lsmeans
##  species lsmean    SE   df lower.CL upper.CL
##  bluti    3.946 0.429 28.7     3.07     4.82
##  greti    0.511 0.763 43.5    -1.03     2.05
##  marti    1.013 1.195 31.4    -1.42     3.45
## 
## Results are averaged over the levels of: Sex, age 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate    SE   df t.ratio p.value
##  bluti - greti    3.434 0.716 48.7   4.795  <.0001
##  bluti - marti    2.932 1.427 29.7   2.055  0.1167
##  greti - marti   -0.502 1.566 32.7  -0.320  0.9451
## 
## Results are averaged over the levels of: Sex, age 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

D. Best model selection

automatic model averaging - with headwind

summary(model.avg(dredge(lmer(v1 ~ species  + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(v1 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(v1 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(v1 ~ species + Sex + age + headwind + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = v1 ~ <9 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##        df  logLik    AICc delta weight
## 4       5 -875.65 1761.34  0.00   0.55
## (Null)  3 -878.91 1763.84  2.50   0.16
## 24      6 -875.91 1763.87  2.53   0.16
## 14      6 -877.12 1766.30  4.96   0.05
## 2       4 -879.64 1767.31  5.96   0.03
## 1       4 -879.73 1767.50  6.15   0.03
## 34      7 -877.36 1768.80  7.46   0.01
## 124     7 -877.41 1768.89  7.55   0.01
## 12      5 -880.45 1770.95  9.60   0.00
## 
## Term codes: 
##      age headwind      Sex  species 
##        1        2        3        4 
## 
## Model-averaged coefficients:  
## (full average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.531459   0.188803    0.188894  13.401   <2e-16 ***
## speciesgreti  0.252291   0.155358    0.155394   1.624    0.104    
## speciesmarti  0.113614   0.128321    0.128404   0.885    0.376    
## headwind     -0.006053   0.013160    0.013162   0.460    0.646    
## age          -0.005385   0.031687    0.031703   0.170    0.865    
## SexF          0.001259   0.022706    0.022721   0.055    0.956    
## SexM          0.002271   0.027944    0.027956   0.081    0.935    
##  
## (conditional average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.53146    0.18880     0.18889  13.401  < 2e-16 ***
## speciesgreti  0.32226    0.09101     0.09109   3.538 0.000403 ***
## speciesmarti  0.14512    0.12830     0.12840   1.130 0.258395    
## headwind     -0.02998    0.01185     0.01186   2.528 0.011466 *  
## age          -0.06034    0.08908     0.08914   0.677 0.498505    
## SexF          0.09467    0.17298     0.17313   0.547 0.584480    
## SexM          0.17075    0.17305     0.17320   0.986 0.324199    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(v2 ~ species  + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(v2 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(v2 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(v2 ~ species + Sex + age + headwind + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = v2 ~ <3 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##    df  logLik    AICc delta weight
## 2   4 -623.79 1255.61  0.00   0.92
## 12  5 -625.32 1260.68  5.07   0.07
## 23  6 -626.57 1265.20  9.59   0.01
## 
## Term codes: 
##      age headwind  species 
##        1        2        3 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.775e+00  1.258e-01   1.258e-01  22.048  < 2e-16 ***
## headwind     -4.589e-02  9.912e-03   9.920e-03   4.626  3.7e-06 ***
## age          -3.349e-03  2.222e-02   2.224e-02   0.151    0.880    
## speciesgreti  2.414e-05  7.129e-03   7.135e-03   0.003    0.997    
## speciesmarti -3.536e-04  1.075e-02   1.075e-02   0.033    0.974    
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.774741   0.125768    0.125848  22.048  < 2e-16 ***
## headwind     -0.045892   0.009912    0.009920   4.626  3.7e-06 ***
## age          -0.046010   0.069432    0.069491   0.662    0.508    
## speciesgreti  0.003173   0.081676    0.081746   0.039    0.969    
## speciesmarti -0.046479   0.114164    0.114261   0.407    0.684    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(avg.v ~ species  + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(avg.v ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(avg.v ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(avg.v ~ species + Sex + age + 
##     headwind + (1 | Ring), data = nat.var, na.action = "na.fail")), 
##     subset = delta < 10)
## 
## Component model call: 
## lmer(formula = avg.v ~ <5 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##        df  logLik    AICc delta weight
## 2       4 -554.97 1117.97  0.00   0.82
## 12      5 -556.04 1122.13  4.16   0.10
## 24      6 -555.73 1123.52  5.54   0.05
## (Null)  3 -559.79 1125.59  7.62   0.02
## 23      6 -557.93 1127.93  9.96   0.01
## 
## Term codes: 
##      age headwind      Sex  species 
##        1        2        3        4 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.639e+00  1.811e-01   1.811e-01  14.569  < 2e-16 ***
## headwind     -3.856e-02  1.073e-02   1.074e-02   3.591 0.000329 ***
## age          -8.302e-03  3.299e-02   3.301e-02   0.252 0.801402    
## speciesgreti  8.398e-03  4.021e-02   4.022e-02   0.209 0.834589    
## speciesmarti  2.826e-03  2.765e-02   2.767e-02   0.102 0.918639    
## SexF          5.814e-05  7.548e-03   7.554e-03   0.008 0.993859    
## SexM          3.090e-04  8.433e-03   8.438e-03   0.037 0.970789    
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.639152   0.181071    0.181144  14.569  < 2e-16 ***
## headwind     -0.039275   0.009445    0.009453   4.155 3.26e-05 ***
## age          -0.080695   0.068838    0.068897   1.171    0.241    
## speciesgreti  0.163409   0.078297    0.078363   2.085    0.037 *  
## speciesmarti  0.054993   0.109572    0.109665   0.501    0.616    
## SexF          0.010274   0.099809    0.099893   0.103    0.918    
## SexM          0.054606   0.097990    0.098073   0.557    0.578    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(acc ~ species  + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(acc ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(acc ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(acc ~ species + Sex + age + headwind + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = acc ~ <8 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##      df   logLik    AICc delta weight
## 34    7 -4141.33 8296.74  0.00   0.18
## 4     5 -4143.35 8296.75  0.01   0.18
## 14    6 -4142.43 8296.92  0.18   0.17
## 134   8 -4140.58 8297.26  0.52   0.14
## 234   8 -4140.97 8298.04  1.30   0.09
## 24    6 -4143.08 8298.21  1.47   0.09
## 124   7 -4142.17 8298.42  1.68   0.08
## 1234  9 -4140.25 8298.62  1.88   0.07
## 
## Term codes: 
##      age headwind      Sex  species 
##        1        2        3        4 
## 
## Model-averaged coefficients:  
## (full average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   5.28324    2.77238     2.77400   1.905   0.0568 .  
## SexF         -0.27265    0.88554     0.88621   0.308   0.7583    
## SexM         -0.49849    0.98070     0.98130   0.508   0.6115    
## speciesgreti -3.40807    0.70124     0.70182   4.856  1.2e-06 ***
## speciesmarti -2.41490    1.23447     1.23543   1.955   0.0506 .  
## age          -0.25223    0.48361     0.48389   0.521   0.6022    
## headwind     -0.06771    0.11703     0.11706   0.578   0.5630    
##  
## (conditional average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)    5.2832     2.7724      2.7740   1.905   0.0568 .  
## SexF          -0.5594     1.2036      1.2046   0.464   0.6424    
## SexM          -1.0228     1.1988      1.1998   0.852   0.3940    
## speciesgreti  -3.4081     0.7012      0.7018   4.856  1.2e-06 ***
## speciesmarti  -2.4149     1.2345      1.2354   1.955   0.0506 .  
## age           -0.5537     0.5886      0.5891   0.940   0.3473    
## headwind      -0.2044     0.1158      0.1159   1.764   0.0778 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

automatic model averaging - with air.temp

summary(model.avg(dredge(lmer(v1 ~ species  + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(v1 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(v1 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(v1 ~ species + Sex + age + air.temp + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = v1 ~ <6 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##        df  logLik    AICc delta weight
## 4       5 -875.65 1761.34  0.00   0.69
## (Null)  3 -878.91 1763.84  2.50   0.20
## 14      6 -877.12 1766.30  4.96   0.06
## 1       4 -879.73 1767.50  6.15   0.03
## 34      7 -877.36 1768.80  7.46   0.02
## 24      6 -879.18 1770.42  9.07   0.01
## 
## Term codes: 
##      age air.temp      Sex  species 
##        1        2        3        4 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.553e+00  1.861e-01   1.862e-01  13.711   <2e-16 ***
## speciesgreti  2.471e-01  1.569e-01   1.570e-01   1.574    0.115    
## speciesmarti  1.117e-01  1.292e-01   1.293e-01   0.864    0.388    
## age          -5.519e-03  3.208e-02   3.210e-02   0.172    0.863    
## SexF          1.566e-03  2.531e-02   2.533e-02   0.062    0.951    
## SexM          2.824e-03  3.114e-02   3.115e-02   0.091    0.928    
## air.temp     -5.167e-05  9.139e-04   9.143e-04   0.057    0.955    
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.552511   0.186079    0.186172  13.711  < 2e-16 ***
## speciesgreti  0.320657   0.091480    0.091557   3.502 0.000461 ***
## speciesmarti  0.145018   0.129793    0.129901   1.116 0.264260    
## age          -0.061645   0.089636    0.089703   0.687 0.491948    
## SexF          0.094674   0.172979    0.173125   0.547 0.584480    
## SexM          0.170751   0.173052    0.173199   0.986 0.324199    
## air.temp     -0.007001   0.008032    0.008039   0.871 0.383801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(v2 ~ species  + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10)) #null AICc higher by 11 delta AIC

## Warning in dredge(lmer(v2 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(v2 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(v2 ~ species + Sex + age + air.temp + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = v2 ~ <3 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##        df  logLik    AICc delta weight
## (Null)  3 -630.73 1267.48  0.00   0.92
## 1       4 -632.26 1272.55  5.07   0.07
## 2       5 -633.45 1276.95  9.47   0.01
## 
## Term codes: 
##     age species 
##       1       2 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.918e+00  1.230e-01   1.230e-01  23.715   <2e-16 ***
## age          -3.223e-03  2.239e-02   2.240e-02   0.144    0.886    
## speciesgreti -5.486e-05  7.511e-03   7.517e-03   0.007    0.994    
## speciesmarti -4.180e-04  1.147e-02   1.148e-02   0.036    0.971    
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.917751   0.122954    0.123032  23.715   <2e-16 ***
## age          -0.044292   0.071203    0.071263   0.622    0.534    
## speciesgreti -0.006806   0.083375    0.083446   0.082    0.935    
## speciesmarti -0.051853   0.116846    0.116945   0.443    0.657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(avg.v ~ species  + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10)) #null AICc higher by 7.62 delta AIC

## Warning in dredge(lmer(avg.v ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(avg.v ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(avg.v ~ species + Sex + age + 
##     air.temp + (1 | Ring), data = nat.var, na.action = "na.fail")), 
##     subset = delta < 10)
## 
## Component model call: 
## lmer(formula = avg.v ~ <4 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##        df  logLik    AICc delta weight
## (Null)  3 -559.79 1125.59  0.00   0.85
## 1       4 -560.88 1129.79  4.20   0.10
## 3       5 -560.83 1131.69  6.10   0.04
## 2       5 -562.63 1135.30  9.71   0.01
## 
## Term codes: 
##     age     Sex species 
##       1       2       3 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.760e+00  1.809e-01   1.810e-01  15.256   <2e-16 ***
## age          -8.309e-03  3.346e-02   3.347e-02   0.248    0.804    
## speciesgreti  6.210e-03  3.440e-02   3.440e-02   0.181    0.857    
## speciesmarti  2.033e-03  2.486e-02   2.487e-02   0.082    0.935    
## SexF          4.366e-05  8.383e-03   8.390e-03   0.005    0.996    
## SexM          3.882e-04  9.488e-03   9.494e-03   0.041    0.967    
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.760498   0.180876    0.180950  15.256   <2e-16 ***
## age          -0.079872   0.071031    0.071091   1.124   0.2612    
## speciesgreti  0.154424   0.080822    0.080890   1.909   0.0563 .  
## speciesmarti  0.050550   0.113628    0.113724   0.444   0.6567    
## SexF          0.006592   0.102790    0.102877   0.064   0.9489    
## SexM          0.058606   0.100899    0.100985   0.580   0.5617    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.avg(dredge(lmer(acc ~ species  + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta <  10))

## Warning in dredge(lmer(acc ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## Warning in dredge(lmer(acc ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML

## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = dredge(lmer(acc ~ species + Sex + age + air.temp + 
##     (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta < 
##     10)
## 
## Component model call: 
## lmer(formula = acc ~ <8 unique rhs>, data = nat.var, na.action = 
##      na.fail)
## 
## Component models: 
##      df   logLik    AICc delta weight
## 34    7 -4141.33 8296.74  0.00   0.25
## 4     5 -4143.35 8296.75  0.01   0.24
## 14    6 -4142.43 8296.92  0.18   0.22
## 134   8 -4140.58 8297.26  0.52   0.19
## 234   8 -4142.53 8301.17  4.43   0.03
## 24    6 -4144.59 8301.24  4.50   0.03
## 124   7 -4143.66 8301.40  4.66   0.02
## 1234  9 -4141.78 8301.69  4.95   0.02
## 
## Term codes: 
##      age air.temp      Sex  species 
##        1        2        3        4 
## 
## Model-averaged coefficients:  
## (full average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   5.445073   2.765768    2.767390   1.968   0.0491 *  
## SexF         -0.274225   0.877139    0.877804   0.312   0.7547    
## SexM         -0.489946   0.970242    0.970840   0.505   0.6138    
## speciesgreti -3.430691   0.697086    0.697663   4.917    9e-07 ***
## speciesmarti -2.421435   1.225460    1.226405   1.974   0.0483 *  
## age          -0.256582   0.484228    0.484500   0.530   0.5964    
## air.temp      0.006876   0.032196    0.032211   0.213   0.8310    
##  
## (conditional average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   5.44507    2.76577     2.76739   1.968   0.0491 *  
## SexF         -0.56945    1.19564     1.19665   0.476   0.6342    
## SexM         -1.01741    1.19087     1.19188   0.854   0.3933    
## speciesgreti -3.43069    0.69709     0.69766   4.917    9e-07 ***
## speciesmarti -2.42144    1.22546     1.22641   1.974   0.0483 *  
## age          -0.56042    0.58469     0.58518   0.958   0.3382    
## air.temp      0.07084    0.07841     0.07848   0.903   0.3667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

final models

lmer.v1.fin <- lmer(v1 ~ species + (1|Ring), data = nat.var)
lmer.v2.fin <- lmer(v2 ~ headwind + (1|Ring), data = nat.var)
lmer.avg.v.fin <- lmer(avg.v ~ headwind + (1|Ring), data = nat.var)
lmer.acc.fin <- lmer(acc ~ species  + Sex + (1|Ring), data = nat.var)

#without outliers
lmer.v1.fin.out <- lmer(v1 ~ species + (1|Ring), data = nat.var.out)
lmer.v2.fin.out <- lmer(v2 ~ headwind + (1|Ring), data = nat.var.out)
lmer.avg.v.fin.out <- lmer(avg.v ~ headwind + (1|Ring), data = nat.var.out)
lmer.acc.fin.out <- lmer(acc ~ species  + Sex + (1|Ring), data = nat.var.out)

E. checking best model assumptions

Acceleration is skipped as the best model was the full model, acceleration only had one outlier, and the full model assumptions were well validated.

normality tests

# Computes the default residuals (raw conditional)
raw_cond.v1.fin <- compute_redres(lmer.v1.fin)
raw_cond.v1.fin.out <- compute_redres(lmer.v1.fin.out)

raw_cond.v2.fin <- compute_redres(lmer.v2.fin)
raw_cond.v2.fin.out <- compute_redres(lmer.v2.fin.out)

raw_cond.avg.v.fin <- compute_redres(lmer.avg.v.fin)
raw_cond.avg.v.fin.out <- compute_redres(lmer.avg.v.fin.out)

#however checking the one without the outliers
raw_cond.acc.fin <- compute_redres(lmer.acc.fin)
raw_cond.acc.fin.out <- compute_redres(lmer.acc.fin.out)

initial velocity

shapiro.test(raw_cond.v1.fin)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v1.fin
## W = 0.96949, p-value < 2.2e-16

shapiro.test(raw_cond.v1.fin.out)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v1.fin.out
## W = 0.98207, p-value = 2.746e-12

end velocity

shapiro.test(raw_cond.v2.fin)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v2.fin
## W = 0.9004, p-value < 2.2e-16

shapiro.test(raw_cond.v2.fin.out)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.v2.fin.out
## W = 0.96104, p-value < 2.2e-16

average velocity

shapiro.test(raw_cond.avg.v.fin)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.avg.v.fin
## W = 0.90085, p-value < 2.2e-16

shapiro.test(raw_cond.avg.v.fin.out)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.avg.v.fin.out
## W = 0.94798, p-value < 2.2e-16

end velocity

shapiro.test(raw_cond.acc.fin)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.acc.fin
## W = 0.9936, p-value = 7.438e-06

shapiro.test(raw_cond.acc.fin.out)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.acc.fin.out
## W = 0.99372, p-value = 1.043e-05

acceleration

shapiro.test(raw_cond.acc.fin)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.acc.fin
## W = 0.9936, p-value = 7.438e-06

shapiro.test(raw_cond.acc.fin.out)

## 
##  Shapiro-Wilk normality test
## 
## data:  raw_cond.acc.fin.out
## W = 0.99372, p-value = 1.043e-05

checking assumptions - no outliers vs. outliers

initial velocity

hist(raw_cond.v1.fin) #normality in histogram

hist(raw_cond.v1.fin.out)

plot_redres(lmer.v1.fin) #heteroskedasticity

plot_redres(lmer.v1.fin.out)

plot_resqq(lmer.v1.fin) #normality of residuals

plot_resqq(lmer.v1.fin.out)

plot_ranef(lmer.v1.fin) #normality of random effects residuals

plot_ranef(lmer.v1.fin.out)

end velocity

hist(raw_cond.v2.fin) #normality in histogram

hist(raw_cond.v2.fin.out)

plot_redres(lmer.v2.fin) #heteroskedasticity

plot_redres(lmer.v2.fin.out)

plot_resqq(lmer.v2.fin) #normality of residuals

plot_resqq(lmer.v2.fin.out)

plot_ranef(lmer.v2.fin) #normality of random effects residuals

plot_ranef(lmer.v2.fin.out)

average velocity

hist(raw_cond.avg.v.fin) #normality in histogram

hist(raw_cond.avg.v.fin.out)

plot_redres(lmer.avg.v.fin) #heteroskedasticity

plot_redres(lmer.avg.v.fin.out)

plot_resqq(lmer.avg.v.fin) #normality of residuals

plot_resqq(lmer.avg.v.fin.out)

plot_ranef(lmer.avg.v.fin) #normality of random effects residuals

plot_ranef(lmer.avg.v.fin.out)

F. Best model summary

initial velocity

#contrasts = list(Ring = contr.sum, cgt_day = contr.sum))
Anova(lmer.v1.fin, type = 3, test.statistic = "F") #species

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v1
##                     F Df Df.res    Pr(>F)    
## (Intercept) 3494.6627  1 39.858 < 2.2e-16 ***
## species        6.3849  2 45.683  0.003588 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lmer.v1.fin.out, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v1
##                     F Df Df.res    Pr(>F)    
## (Intercept) 3720.5816  1 39.989 < 2.2e-16 ***
## species        6.1535  2 45.552  0.004309 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.v1.fin)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v1 ~ species + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1751.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7220 -0.5874  0.0470  0.6827  3.5530 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.04864  0.2205  
##  Residual             0.18537  0.4305  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   2.50680    0.04235 38.68447  59.190  < 2e-16 ***
## speciesgreti  0.32179    0.09112 52.30977   3.531 0.000872 ***
## speciesmarti  0.14344    0.12596 38.05612   1.139 0.261922    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr
## speciesgret -0.465       
## speciesmart -0.336  0.156

tab_model(lmer.v1.fin)

	v1
Predictors	Estimates	CI	p
(Intercept)	2.51	2.42 – 2.59	<0.001
species [greti]	0.32	0.14 – 0.50	<0.001
species [marti]	0.14	-0.10 – 0.39	0.255
Random Effects
σ²	0.19
τ₀₀ _Ring	0.05
ICC	0.21
N _Ring	49
Observations	1434
Marginal R² / Conditional R²	0.042 / 0.241

lsmeans(lmer.v1.fin, pairwise~species)

## $lsmeans
##  species lsmean     SE   df lower.CL upper.CL
##  bluti     2.51 0.0424 39.9     2.42     2.59
##  greti     2.83 0.0809 59.0     2.67     2.99
##  marti     2.65 0.1187 39.1     2.41     2.89
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate     SE   df t.ratio p.value
##  bluti - greti   -0.322 0.0913 53.9  -3.524  0.0025
##  bluti - marti   -0.143 0.1260 39.2  -1.138  0.4967
##  greti - marti    0.178 0.1436 44.2   1.242  0.4352
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

lsmeans.v1 <- lsmeans(lmer.v1.fin, pairwise~species)

end velocity

Anova(lmer.v2.fin, type = 3, test.statistic = "F") #species and headwind

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v2
##                    F Df  Df.res    Pr(>F)    
## (Intercept) 3857.990  1  160.93 < 2.2e-16 ***
## headwind      21.387  1 1431.16 4.091e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lmer.v2.fin.out, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: v2
##                    F Df  Df.res    Pr(>F)    
## (Intercept) 4472.346  1  142.06 < 2.2e-16 ***
## headwind      21.269  1 1413.59 4.351e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.v2.fin)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v2 ~ headwind + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1247.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8745 -0.4077  0.0902  0.5318  3.5448 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.03895  0.1973  
##  Residual             0.12968  0.3601  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  2.757e+00  4.433e-02  1.571e+02   62.19  < 2e-16 ***
## headwind    -4.589e-02  9.912e-03  1.431e+03   -4.63 3.99e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## headwind 0.700

average velocity

Anova(lmer.avg.v.fin, type = 3, test.statistic = "F") #species and headwind

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: avg.v
##                    F Df  Df.res    Pr(>F)    
## (Intercept) 3587.828  1  149.85 < 2.2e-16 ***
## headwind      17.223  1 1429.98  3.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lmer.avg.v.fin.out, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: avg.v
##                    F Df  Df.res    Pr(>F)    
## (Intercept) 4221.690  1  135.54 < 2.2e-16 ***
## headwind      13.393  1 1412.39 0.0002619 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.avg.v.fin)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.v ~ headwind + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 1109.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6928 -0.4278  0.1154  0.6222  3.7496 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept) 0.03948  0.1987  
##  Residual             0.11744  0.3427  
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  2.595e+00  4.328e-02  1.428e+02  59.969  < 2e-16 ***
## headwind    -3.924e-02  9.444e-03  1.430e+03  -4.155 3.45e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## headwind 0.683

acceleration

Anova(lmer.acc.fin, type = 3, test.statistic = "F") #species and Sex

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: acc
##                   F Df Df.res    Pr(>F)    
## (Intercept) 16.0819  1 28.151 0.0004057 ***
## species     12.7559  2 39.110 5.436e-05 ***
## Sex          0.6363  2 33.651 0.5354960    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lmer.acc.fin.out, type = 3, test.statistic = "F")

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: acc
##                  F Df Df.res    Pr(>F)    
## (Intercept) 15.522  1 28.253 0.0004878 ***
## species     13.214  2 39.044 4.141e-05 ***
## Sex          0.493  2 33.764 0.6150967    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lmer.acc.fin)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: acc ~ species + Sex + (1 | Ring)
##    Data: nat.var
## 
## REML criterion at convergence: 8282.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3766 -0.6433  0.0123  0.6439  3.4045 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Ring     (Intercept)  2.162   1.470   
##  Residual             18.181   4.264   
## Number of obs: 1434, groups:  Ring, 49
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    4.4944     1.1206 27.0775   4.011 0.000429 ***
## speciesgreti  -3.2851     0.6995 50.4970  -4.696 2.07e-05 ***
## speciesmarti  -2.6814     1.4065 29.5907  -1.906 0.066347 .  
## SexF          -0.5125     1.1964 27.9608  -0.428 0.671682    
## SexM          -1.0269     1.1977 28.0081  -0.857 0.398506    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) spcsgr spcsmr SexF  
## speciesgret  0.000                     
## speciesmart -0.797  0.000              
## SexF        -0.937 -0.089  0.746       
## SexM        -0.936 -0.149  0.745  0.891

tab_model(lmer.acc.fin)

	acc
Predictors	Estimates	CI	p
(Intercept)	4.49	2.30 – 6.69	<0.001
species [greti]	-3.29	-4.66 – -1.91	<0.001
species [marti]	-2.68	-5.44 – 0.08	0.057
Sex [F]	-0.51	-2.86 – 1.83	0.668
Sex [M]	-1.03	-3.38 – 1.32	0.391
Random Effects
σ²	18.18
τ₀₀ _Ring	2.16
ICC	0.11
N _Ring	49
Observations	1434
Marginal R² / Conditional R²	0.062 / 0.161

lsmeans(lmer.acc.fin, pairwise~species)

## $lsmeans
##  species lsmean    SE   df lower.CL upper.CL
##  bluti    3.981 0.429 29.6    3.105     4.86
##  greti    0.696 0.737 46.9   -0.787     2.18
##  marti    1.300 1.151 32.4   -1.044     3.64
## 
## Results are averaged over the levels of: Sex 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate    SE   df t.ratio p.value
##  bluti - greti    3.285 0.702 52.5   4.679  0.0001
##  bluti - marti    2.681 1.407 30.8   1.906  0.1541
##  greti - marti   -0.604 1.572 33.8  -0.384  0.9221
## 
## Results are averaged over the levels of: Sex 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

lsmeans.acc <- lsmeans(lmer.acc.fin, pairwise~species)

plotting species difference

plot.v1.means <- data.frame(lsmeans.v1$lsmeans)


plot.acc.means <- data.frame(lsmeans.acc$lsmeans)

plot.v1.lsmeans <- 
  ggplot() +  
  geom_point(data = plot.v1.means, aes(x = species, y = lsmean, color = species), alpha = 1) +
  geom_errorbar(data = plot.v1.means, aes(x = species, ymin = lower.CL, ymax = upper.CL, color = species, width = 0.5),alpha = 1) +
  scale_color_manual(values = c("bluti" ="black", "greti" = "black", "marti" = "black")) +
  scale_x_discrete(labels = c("Blue Tit","Great Tit","Marsh Tit")) +
  ylab(expression(initial~velocity~(m/s))) +
  ylim(2.25,3.25)+
  theme(axis.title.x=element_blank()) + #getting rid of x-axis title
  theme(legend.position = "none")  #getting rid of legend

plot.acc.lsmeans <- 
  ggplot() +  
  geom_point(data = plot.acc.means, aes(x = species, y = lsmean, color = species), alpha = 1) +
  geom_errorbar(data = plot.acc.means, aes(x = species, ymin = lower.CL, ymax = upper.CL, color = species, width = 0.5),alpha = 1) +
  scale_color_manual(values = c("bluti" ="black", "greti" = "black", "marti" = "black")) +
  scale_x_discrete(labels = c("Blue Tit","Great Tit","Marsh Tit")) +
  ylab(expression(`in`-flight~acceleration~(m/s^2))) +
  ylim(-2,6) +
  theme(axis.title.x=element_blank()) + #getting rid of x-axis title
  theme(legend.position = "none") #getting rid of legend

grid.arrange(plot.v1.lsmeans, plot.acc.lsmeans, nrow = 1)

Natural variation of flight

Kyu Min Huh

2024-06-18

Introduction

I. Data

libraries

dataset

adding day and hour and climate

data synopsis

raw continuous data

natural variation data

For each species, M vs. F, young vs. adult, then between species

no.of individuals by species, age, sex in natural variation

t-tests on visiting rate

headwind, air temperature, humidity

duration difference

II. Graphs

A. Structural variation

1. individual variation

2. species variation

3. sex x species variation

4. age variation

B. climate factors

5.1 headwind

5.2 air temperature

5.3 humidity

C. diurnal variation

III. Linear Mixed Model analysis

A. Model selection

1. checking correlations

2. Proposed Mixed Linear Models.

with headwind

with air temperature

3. preparing dataset without outliers

outliers

without outliers

B. Linear mixed model assumptions

normality tests

initial velocity

end velocity

average velocity

end velocity

checking assumption graphs

initial velocity

end velocity

average velocity

acceleration

C. Model summary

initial velocity

end velocity

average velocity

acceleration

D. Best model selection

automatic model averaging - with headwind

automatic model averaging - with air.temp

final models

E. checking best model assumptions

normality tests

initial velocity

end velocity

average velocity

end velocity

acceleration

checking assumptions - no outliers vs. outliers

initial velocity

end velocity

average velocity

F. Best model summary

initial velocity

end velocity

average velocity

acceleration

plotting species difference